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The carrier density distributions in few-layer-graphene systems grown on the carbon face of silicon 
carbide can be altered by the presence of a scanning tunneling microscope (STM) tip used to probe 
top-layer electronic properties, and by a perpendicular magnetic field which induces well-defined 
Landau levels. Hartree approximation calculations in the perpendicular field case show that charge 
tends to rearrange between the layers so that the filling factors of most layers are pinned at integer 
values. We use our analysis to provide insight into the role of buried layers in recent few-layer- 
graphene STM studies and discuss the limitations of our model. 

PACS numbers: 73.22.Pr, 68.37.Ef, 71.70.Di 



I. INTRODUCTION 

Progress in the preparation and isolation of highly 
ordered graphene sheets over the past few yearsi 2 - has 
led to an explosion of interest in the properties of these 
two-dimensional electron systems which are remarkably 
simple, yet rich in interesting mechanical and electronic 
properties. One type of graphene system 2 - that is po- 
tentially suitable for applications is prepared by thermal 
decomposition of silicon carbide (SiC). The unique fea- 
ture of these epitaxial graphene systems is that they tend 
to grow not as single layers but as few layer graphene 
(FLG) systems. The layers tend to be electrically isolated 
to a reasonable degree^— because of partially controlled 
relative rotations^ FLG systems on SiC can be grown 
as large area films that are extremely highly ordered, at 
least locally, and doped by charge transfer from a carbon 
buffer layer, which is a nongraphitic carbon layer between 
the SiC and the graphene layers. 

This paper addresses the distribution of charge carriers 
across the FLG system. Recent measurements of Landau 
level spectr a 10 ' 11 and angle-resolved photoemissiooi 2 - for 
FLG systems grown on the carbon face of SiC show char- 
acteristics of decoupled monolayer graphene rather than 
coupled graphene multilayers. This behavior is likely due 
to the relative rotations between the layers. Our work is 
motivated in part by an interest in understanding the 
confusing^ transport properties of these systems, which 
must be strongly dependent on carrier charge distribu- 
tion across the weakly coupled layers. Our immediate 
motivation, however, is provided by recen t 11 ' 14 scanning 
tunneling microscopy (STM) Landau level spectroscopy 
studies of FLG in the presence of an external magnetic 
field. Traces of the Landau level positions can be ex- 
tracted from such spectra as is illustrated in Fig. [TJ Al- 
though STM directly probes electronic properties in the 
top layer of a FLG system, there is evidence that top 
layer properties can be altered, sometimes qualitatively, 
by correlations with electrons in submerged layers. 

If the density in the top layer were fixed, the Fermi 
level would be pinned to one of the Landau level energies 
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FIG. 1: (Color online) Landau level peak positions of epi- 
taxial graphene on C-face SiC as a function of magnetic field 
obtained from the STM measurements in Ref. The posi- 
tion of each Landau level is averaged over spin and valley split 
dl/dV peaks when these can be separately resolved.— dl/dV 
peak positions at finite tip-sample bias can be influenced by 
tip-sample band bending and by redistributions of charge in 
the FLG system as explained in detail below. 

except at the discrete field strengths which yield integer 
filling factors when the Fermi level is in between Landau 
levels. In practice, experiment shows the opposite behav- 
ior. The Landau levels tend to be pinned away from the 
Fermi energy, an effect that is particularly striking in the 
field range between 8 T and 10 T. At higher fields, the 
Landau levels split through valley and spin splitting so 
that even in the field range near 12 T, the split Landau 
levels avoid the Fermi energy. 

The most dramatic effect seen in these experiments 
is splitting within spin- and valley-split iV=l peaks in 
the density of states (DOS) as they pass through the 
Fermi leveli 1 ^ This peculiar, fractionally filled Landau 
level gives evidence for a correlated-electron state that is 
stable when the N = 1 Landau level of the top layer is 
half-filled. 

While the precise nature of this state remains myste- 
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rious, its formation might depend only on correlations 
among top layer electrons; however, if one of the sub- 
merged layers is also partially filled under the same tip- 
biasing and field conditions, then this fractionally filled 
Landau level could depend essentially on correlations be- 
tween electrons in different layers. Since half-filling does 
not favor the formation of especially stable states in an 
isolated layer, the latter possibility appears likely. In 
the strong-magnetic-field quantum Hall regime with fully 
formed Landau levels only weakly broadened by disor- 
der, correlations are strongest when Landau levels are 
partially filled. One goal of the model developed in this 
paper attempts to provide a basis for estimating which 
layers contain partially filled Landau levels as the mag- 
netic field strength varies. 

Our paper is organized as follows. In Sec. II, we ex- 
plain our model for carrier distribution in a few layer 
graphene system in which the buffer layer acts as a 
reservoir for carriers. We assume that the growth- 
dependent buffer layer properties determine the position 
of the Fermi energy relative to the Dirac point of the 
bottom graphene layer. Electron-electron interactions 
are included only at the Hartree level. At zero mag- 
netic field carriers reside mainly in the layers closest to 
the buffer and the density-of-states in the top layer is 
small. When a perpendicular magnetic field is applied, 
the Fermi level tends to be pinned near one of the filling 
factors (y = ±2, ±6, ±10, . . .) at which the integer quan- 
tum Hall effect occurs in a graphene layer, implying that 
charge must be transferred between layers as a function 
of field. In Sec. Ill, we discuss how an STM tip can be 
included in such a model. When an STM tip is intro- 
duced to study the electronic properties of the top layer, 
its carrier density tends to be altered with a sign and 
magnitude that is strongly dependent on the tip work 
function. The STM studies in Refs. show that the 

top-layer is n-type for the tip used in those experiments, 
so that carrier densities peak not only near the buffer 
layer but also near the top-layer. In Sec. IV, we use this 
basic theoretical picture to develop a theory of STM Lan- 
dau level spectroscopy in FLG, comparing where possible 
with STM data. We find that as the sample-tip bias and 
the magnetic field are varied, charge tends to rearrange 
to achieve integer filling factors in as many of the FLG 
layers as possible. In Sec. V, we conclude with a brief 
summary and some suggestions for future experimental 
and theoretical work. 



II. FEW-LAYER-GRAPHENE MODEL 

We estimate carrier charge distribution in a FLG sys- 
tem grown on carbon-face SiC substrates using the model 
summarized schematically in Fig. [5J Earlier work con- 
sidered the charge distribution on mono and bilayer 
graphene^ and for multilayer graphene- 7 in the contin- 
uum limit, both in the absence of a magnetic field, our 
main interest. In Fig. [3J the graphene layers are labeled 
by integer numbers starting from label 1 for the layer 




FIG. 2: (Color online) Schematic illustration of a few-layer- 
graphene system without an STM tip. The buffer layer be- 
tween the SiC and the graphene layers was omitted for sim- 
plicity. 

closest to the buffer layer to M for the top layer. (The 
buffer layer between the SiC and the graphene layers, 
which acts as a reservoir for carriers, was omitted for sim- 
plicity.) In equilibrium all layers share the same chemical 
potential \x. The Dirac point in layer i is shifted by its 
local electric potential Ui. The energy spectrum of each 
layer is that of monolayer graphene with the Dirac point 
shifted by Ui. It follows that the charge density in layer 
i, Ui, satisfies 

^ _ sgn(/x - in) ( fi - Ui \ 2 

7T \ Hv J 

The potential energy ttj is in turn evaluated from the 
charge densities using the Poisson equation which implies 
that the electric field Ei between layer i and layer i + 1 
satisfies 

e(Ei-E i _ 1 )=4w{-e)a i . (2) 

The dielectric constant e in Eq. @ accounts for the polar- 
izability between graphene sheets. Here we choose e = 1; 
we have found that changing the value of e does not qual- 
itatively alter the main results of this paper. 

Since the electric field in the vacuum above the top 
(M-th) layer E vac must vanish in the absence of an STM 
tip, the electric fields between all graphene sheets are 
readily evaluated iteratively given the charge densities 
Ui. Starting from layer 1 and adding a contribution due 
to the electric field between a layer and the layer above 
gives: 

u i+ i = Ui + edEi, (3) 

where d — 0.335 nm is the interlayer separation between 
graphene layers. Our neglect of interaction effects be- 
yond electrostatics is supported in the zero magnetic field 
limit by recent Green's function screened Coulomb (GW) 
many-body calculations^ by Profumo et al. As we dis- 
cuss below, exchange and correlation effects are likely to 
be more important in large magnetic fields. 

We model the role of the buffer layer by assuming that 
its equilibration with the bottom graphene layer fixes the 
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FIG. 3: (Color online) Charge density vs. layer index for a 
6-layer FLG system at B = T and B = 10 T. Here chemical 
potential fi — 360 meV and temperature T — 30 K were used. 



value of fj, — ui . It is known that the carrier density of the 
graphene system is sensitive to the microstructure of the 
disordered buffer layer, and hence to FLG growth con- 
ditions. For a given sample, some carriers remain after 
the bonding between the buffer layer and the SiC sub- 
strate is established. Energy in the system is lowered as 
these electrons are transferred to the 7r-bands of the first 
graphene layer. u± is determined by a balance between 
the chemical driving force for the transfer and the band 
and electrostatic energy cost of adding electrons to the 
graphene. In assuming that u\ is independent of field, 
as we do below, we are taking advantage of the fact that 
the Landau-level separation at the Fermi energy of the 
first graphene sheet is small compared to fi — u±. In 
modeling STM data on a particular sample, imperfect 
knowledge of the most appropriate value for fi — u\ is 
an important source of uncertainty that limits predictive 
power. For the calculations described below we choose 
fi—u\ = 360 meV, an estimate that is motivated by spec- 
troscopic measurements^ in multilayer graphene grown 
on the C-face of the SiC substrate. In the rest of this 
paper we choose our zero of energy so that u\ = 0. 

To explain how band and electrostatic energies com- 
bine to determine carrier distributions, we first consider 
double-layer graphene with a chemical potential fi > 0. 
From the Poisson equation with E waE = 0, the electric 
field between layers, E, satisfies eE — \-Kea 2 . The po- 
tential energy of the top layer is u 2 = eEd. It follows 
that 
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<?2 = 



sgn(fi ~u 2 ) ( ' fi — u 2 



7T \ hv 

sgn(n - u 2 ) fJJ_\ 2 
hv 



where 



and (3 = e(hv) 2 /(8e 2 dfi) is a unitless quantity which con- 
trols the energy balance of charge moving between layers. 
Here v is the 7T-band velocity at the Dirac point which 
is proportional to the intralayer hopping energy 70. For 
fi = 360 meV, e = 1, and 70 = 3 eV, [3 ~ 0.29 and 
/(/?) fa 0.28; thus the top layer has 28% of the bottom 
layer charge. 

For multilayers with more than two layers, the distri- 
bution follows from a simple numerical calculation. The 
layer charge density is calculated by integrating the Lan- 
dau level density of states weighted by the Fermi factor 
for the appropriate chemical potential. Then from the 
resulting layer densities, layer potentials are calculated 
using the Poisson equation in Eq. (jSJ. This process is 
repeated until a self-consistency is reached. 

The charge distribution for a decoupled 6-layer 
graphene stack at B = T and B = 10 T calculated with 
the same parameters used in the double layer graphene, 
is shown in Fig. [31 For B = T, layers above the bot- 
tom layer have in total 32% of the bottom layer charge 
and this ratio is almost independent of the number of 
layers. As the magnetic field is turned on, the charge 
distribution is altered due to the formation of Landau 
levels, particularly in the low-density layers with a Fermi 
level near the Dirac point. Because a Landau level ap- 
pears precisely at the Dirac point in a graphene sheet, a 
magnetic field causes a peak in the density-of-states to 
appear at the same energy at which the density-of-states 
vanishes in the absence of a magnetic field. This feature 
of graphene physics strengthens magnetoelectric effects 
associated with Landau level quantization. 



III. STM TIP MODEL 

Figure 0] shows a schematic illustration of a FLG sys- 
tem with an STM tip. We model an STM tip as an 
additional layer which acts as a top gate electrode. The 
distance between the tip and graphene surface is taken 
as civac = 1 nmi2 

Experimentally, it is found that the graphene work 
function, i.e., the energy to take an electron from the 
Fermi energy to vacuum, depends on the charge on the 
surface layer. However, it is also found that the en- 
ergy to take a graphene electron from the Dirac point 
to vacuum does not change as a function of the charge 
density— We denote the latter as <J? gr . In general, the 
work function of the tip (3>tip) and the graphene layer 
[$ gr — (/1 — um)] are different. This difference in work 
functions, 'f'gr — (fi — um) — ^tip = $ — (a* ~ u m) leads to 
charge transfer between the surfaces when they are elec- 
trically connected and induces an electric field between 
the surface and tip. As seen in the bottom panel of Fig. 0] 
the electric field satisfies 



fi + eV = fi 



tip 



U M 



(6) 



(5) 



as a voltage V is applied between tip and sample. 

In STM spectroscopy a new tunneling transport chan- 
nel opens up, giving rise to a dl/dV peak whenever the 
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FIG. 4: (Color online) Upper panel shows schematic illustra- 
tion of a few-layer-graphene system with an STM tip. The 
buffer layer between the SiC and the graphene layers was 
omitted for simplicity. Lower panel is the energy level dia- 
gram for this system. 



chemical potential of the STM tip is aligned with one 
of the top- layer Landau levels. The experimental dl/dV 
peaks therefore identify the tip-sample bias voltages at 
which the following resonant tunneling conditions are sat- 
isfied: 



fi + eV = u M + e N (B), 



(7) 



where En{B) — sgn(N) y '2\N\ehv 2 B / 'c is the graphene 
sheet Landau level energy. To illustrate the effect of the 
tip, in Fig. [5] we calculate the charge distribution of each 
layer at B = 10 T for $ = +0.4 eV, eV, and -0.4 eV 
when the tip-sample bias V is zero. 

For $ = eV and V = V, the FLG charge distribu- 
tion is identical to the distribution without an STM tip 
at B — 10 T shown in Fig. [3l For non-zero <£, however, 
an electric field between the tip and sample surface is in- 
duced and distorts the layer charge distribution even at 
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FIG. 5: (Color online) Charge density in the layers of a six 
layer stack at B = 10 T for $ = +0.4 eV, eV, and -0.4 eV 
when the tip-sample bias is zero. These curves were obtained 
using /i = 360 meV and T = 30 K. By construction, electrons 
have a positive carrier density, while holes have negative den- 
sity (see text for details). 



zero tip-sample bias. 



IV. FLG LANDAU-LEVEL TUNNELING 
SPECTROSCOPY 

At weak fields, um and E vac are approximately con- 
stant so that the spacing in electronvolts between dl/dV 
peaks matches the energetic separation between top-layer 
Landau levels. The spectroscopy data can therefore be 
used to measure the Dirac velocity parameter which char- 
acterizes the energy scale of the graphene layer's Dirac 
cones. In the strong-field limit, however, the density-of- 
states in each graphene layer is altered, and this in turn 
alters the densities at which equilibria are established be- 
tween adjacent layers. It follows that both um and -E V ac 
depend on field. One goal of our calculations is to esti- 
mate the magnitude and character of this effect. 

To illustrate this effect we first examine the field- 
dependence of the N = Dirac-point Landau level fea- 
ture, plotted in Fig. O Since En=o(B) = 0, the field- 
dependence of this spectral feature is due entirely to the 
field-dependence of charge distributions in the FLG sys- 
tem. The calculations in Fig. [5] were carried out at a 
finite temperature T = 30 K, in part to crudely model 
the Landau-level smearing influence of disorder. At weak 
fields the N = dl/dV peak's position is independent of 
field as expected. The position of these peaks is primarily 
dependent on the model's workfunction parameter $. 

As illustrated in Fig. [6j the position of the experimen- 
tal peak at eV « —135 meV and at B = 5 T is re- 
produced approximately by setting /i = 360 meV, which 
leads to u M ~ 225 meV and $ w 400 meV. The slow 
downward drift in the weak field N — Dirac point 
peak with increasing magnetic field is not reproduced by 
our calculation, and could be due to an increase in the 
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FIG. 6: (Color online) The field dependence of the iV = 
Dirac-point Landau level for different workfunction parame- 
ters <£>.— Here fj, = 360 meV and T = 30 K were used. 

strength of exchange and correlation effects in FLG with 
magnetic field. The strong variations in peak positions 
with field that begin at around 6 T are the quantizing- 
magnetic-field effects on which we will focus in the re- 
mainder of this paper. 

The influence of Landau quantization on dl/dV spec- 
tra is illustrated in more detail in Fig. which shows the 
prediction of the theoretical model for $ = 400 meV and 
H = 360 meV at T = 30 K. 

When the Landau level energies in a particular layer 
are far away from the Fermi energy, the layer filling 
factor Vi — 2nt 2 (Ti, where i = yj hc/eB is pinned at 
one of the full-Landau-level filling factor values: vi — 
±2,±6,±10,---. For a fixed filling factor the carrier den- 
sity in a layer increases with field and its Landau level 
energies therefore increase due to electrostatic repulsion. 
The increase in density must be achieved by charge trans- 
fer from other layers. When a Landau level in a layer is 
close to the Fermi level, the density in that layer will tend 
to decrease as its Landau level empties with increasing 
field. This is the source of charge transferred to other 
layers. This behavior contrasts with that of an isolated 
system with fixed charge density in which integer filling 
factors occur only at isolated field values and successive 
Landau-level energies are pinned to the Fermi level. 

We refer to layers which have partially filled Landau 
levels as active and to layers which have full Landau levels 
as inactive. Since the total filling factor is a smooth func- 
tion of field, at least one layer must be active at generic 
field values. Strong interlayer correlations are likely when 
two or more layers are active. It would be surprising if 
interlayer correlation effects were not important, given 
the relationship between the important length scales in 
the problem. At 10 T, the total graphene layer thickness 
Md « 2 nm is typically much less than the average sep- 
aration of electrons within one layer l/^/al ~ 10 nm as 
well as the magnetic length I — y/hc/eB « 8.1 nm. Such 
correlations could give rise to a state with spontaneous in- 
terlayer coherence^ among other possibilities. As shown 
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FIG. 7: (Color online) Theoretical prediction of top layer 
Landau-level peak positions as a function of magnetic field 
for $ = 400 meV and fi = 360 meV at T = 30 K. 

below, our model calculations provide estimates of the 
field ranges at which two or more layers become active. 

For the parameters of Fig. [7] the model predicts that 
the top layer is active between B = 5 T and B = 7 T. In 
this field range the N = 2 Landau level is pinned to the 
Fermi level and the filling factor varies between v = 10 
and v = 6. 22 The top layer is then briefly inactive before 
becoming active again above 8 T when the N = 1 Lan- 
dau level is pinned to the Fermi level. In inactive field 
ranges, the density in the top layer is proportional to 
magnetic field and the energies of all levels in that layer 
increase. In the active regions, the density tends to de- 
crease and the rate at which energy levels increase with 
field for N > is suppressed by the decrease in density. 
This is only a tendency, however, since the electrostatic 
energy in the top layer depends on the densities at all 
layers. The evolution of the STM spectrum also depends 
on the evolution of charge density in the submerged lay- 
ers that are not directly probed by the STM. Note that 
the distribution of charge among the FLG layers also de- 
pends somewhat on the tip-sample bias voltage. Figure [8] 
shows the electrostatic energy, density, and filling factor 
of each layer as a function of magnetic field at zero bias 
voltage. These results are consistent with the preceding 
discussion. 



V. DISCUSSION AND CONCLUSIONS 

Our theoretical model does not account for exchange 
and correlation effects, which can alter the energy change 
associated with adding electrons to empty (or partially 
filled) Landau levels and the energy change associated 
with removing electrons from full (or partially full) Lan- 
dau levels. Systematic discrepancies between present the- 
ory and experiment likely signal these neglected interac- 
tion effects. These discrepancies include the low field 
variation of the zcroth Landau level energy (Fig. [5]) and 
the pinning of the first Landau level away from the Fermi 
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FIG. 8: (Color online) Electrostatic energy, density and filling 
factor of each layer as a function of magnetic field for $ = 400 
meV, /i = 360 meV and T = 30 K at zero tip-sample bias 
voltage. The filling factor of the bottom layer (not shown) 
exceeds v = 12 over the field range considered. Numbers 
indicate layer numbers, with 6 being the topmost layer. 



energy rather than at it (Fig. [TJ. 

In Fig. IH1 for example, we compare experimental and 
theoretical energy separations between N = and the 
N = 1,2 dl/dV features as a function of magnetic field 
strength. Even at the Hartree level there arc additions to 
the \/B band energy contribution due to changes in elec- 
trostatic energies with tip-sample bias voltages indicating 
ovcrestimation of the electrostatic tip-gating effects. 

We note that the filling-factor dependent features 
in the field-dependence of the Landau level energies 
[Fig. EJc)] are weaker in experiment than in this Hartree 
theory. We believe that these differences mainly reflect 
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FIG. 9: (Color online) Landau level separations between 
N = and the N = 1 (N = 2), denoted as LLi (LL 20 ) 
for theoretical calculations and STM measurements.— For the 
theoretical calculations, "!> = 400 meV and /i = 360 meV at 
T = 30 K were used. 
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FIG. 10: (Color online) Filling factor modulus i^ mod as a func- 
tion of magnetic field for $ = 400 meV and fi = 360 meV 
at T = 30 K, where v™ od = |(z//4 mod 1) - 1/2|. Note 
that vr od = for v t =±2, ±6, ±10, ■ ■ ■ , while v™ oi = 1/2 for 
Vi—Q, ±4. ± 8, ■ • ■ . The points were evaluated at zero bias 
voltage. 

exchange and correlation energies which mitigate elec- 
trostatic effects. The presence of strong correlation ef- 
fects in this field range is apparent in the experimental 
interaction-induced Landau level splittings which have 
been suppressed in Fig. |H] by averaging over all experi- 
mental features identified with 7V=0 and N=l. 

The approximately linear reduction with field (at weak 
fields) of the bias voltage at which the N = peak is 
observed (see Fig. [5]) is completely absent in theory and 
unexplained at present. 

In addition to effects associated with exchange and cor- 
relation within a layer, we expect that interlayer corre- 
lations play an essential role when more than one layer 
is active, at least when the magnetic field is strong and 
Landau levels are well developed. In Fig. HUl we plot the 
field dependence of the partial filling factor (per layer 
and spin) of the active Landau level: ^™ od = |(f/4 mod 
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1) — 1/2 1, which we refer to as the modular filling factor. 
v mod | g defined so that it vanishes when a Landau level 
is either completely filled or completely empty, namely 
at total filling factors z/j=±2, ±6, ±10, • • • . Note that the 
layer is inactive when v™- od = 0, and most active when 
v™° d = 0.5. Around B = 11 T, we see that the top (6th) 
layer and the 2nd layer have filling factors close to half- 
filled filling factors, v§ — 4 and v 2 = 8, respectively (in 
this field range the intervening layers all have total Lan- 
dau level filling factor v = 2). This is precisely the field 
range in which a gap appears to open in the top layer 
N = 1 tunneling density of states in STM studies^ The 
gap could therefore be due to correlations between N = 1 
electrons in the top layer and N = 2 electrons in layer 
2. The participation of electrons in another layer could 
explain the appearance of a gap at a partial filling factor 
which is not known to support large gaps in a single- 
layer system. One possible state that is consistent with 
experiments is one in which coherence is spontaneously 23 
established between layers 2 and 6. Around B = 5 T and 
below, the top layer and 5th layers become active, but 



due to the smaller field magnitude, strong correlation ef- 
fects are more easily suppressed by disorder. 

The analysis presented in this paper highlights both 
advantages and disadvantages of few layer graphene sys- 
tems for physics studies. Because the electronic degrees 
of freedom in all layers can play an active role, particu- 
larly at strong magnetic fields, the physics is extremely 
rich. On the other hand, the same property makes it 
more challenging to uniquely interpret observations us- 
ing surface physics probes, like STM, which are directly 
sensitive mainly to top layer properties. 
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